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EXPONENTIAL APPROXIMATIONS USING FOURIER SERIES PARTIAL SUMS 

NANA S. BANERJEE* AND JAMES F. GEERt 

Abstract. The problem of accurately reconstructing a piece- wise smooth, 27r~pcriodic function / and 
its first few derivatives, given only a truncated Fourier series representation of /, is studied and solved. The 
reconstruction process is divided into two steps. In the first step, the first 2N + 1 Fourier coefficients of / 
are used to approximate the locations and magnitudes of the discontinuities in / and its first M derivatives. 
This is accomplished by first finding initial estimates of these quantities based on certain properties of Gibbs 
phenomenon, and then refining these estimates by fitting the asymptotic form of the Fourier coefficients to 
the given coefficients using a least-squares approach. It is conjectured that the locations of the singularities 
are approximated to within O (JV”^ -2 ) , and the associated jump of the k th derivative of / is approximated 
to within O (jV _M_1+ *) , as N — ► oo, and the method is robust. These estimates arc then used with a 
class of singular basis functions, which have certain “built-in” singularities, to construct a new sequence of 
approximations to /. Each of these new approximations is the sum of a piecewise smooth function and a 
new Fourier series partial sum. When N is proportional to M, it is shown that these new approximations, 
and their derivatives, converge exponentially in the maximum norm to /, and its corresponding derivatives, 
except in the union of a finite number of small open intervals containing the points of singularity of /. The 
total measure of these intervals decreases exponentially to zero as M — * oo. The technique is illustrated with 
several examples. 

Key words. Fourier series, exponentially accurate approximations, piecewise smooth functions, location 
of singularities 

Subject classification. Applied and Numerical Mathematics 

1. Introduction. Approximate solutions to problems in applied mathematics arc often obtained using 
a finite number of terms in the Fourier series representation of the solution. In practice, this truncation 
procedure may lead to nonuniformly valid approximations. In particular, when the function being approx- 
imated has one or more points of discontinuity, Gibbs phenomena is present, resulting in an “overshoot” 
of the jump in the function at a point of discontinuity, as well as artificial oscillations near such a point. 
The magnitude of the overshoot is not eliminated by increasing the number of terms in the approximation. 
In addition, the oscillations caused by this phenomena typically propagate into regions away from the sin- 
gularity, and, hence, degrade the quality of the partial sum approximation in these regions. It has been 
conjectured, however, that this oscillatory approximation, which may have been obtained by a high-order 
method, such as a spectral method, should contain enough information to enable the reconstruction of the 
proper, non-oscillatory, discontinuous function, by a post-processing filter (see, e.g., Lax [17]). 

In a series of papers, Gottlieb, et. al. [7], [8], [9], [10], [11], have proposed and investigated a way 
of overcoming the Gibbs phenomena. Their technique involves the construction of a new series using the 
Gegenbauer polynomials. For a function / that is analytic on the interval [—1,1], but is not periodic, they 
prove that their technique leads to a series which converges exponentially to / in the maximum norm. Re- 
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cently, Geer [5] introduced and studied a class of approximations {FW,m} to a periodic function / which 
uses the ideas of Pade approximations based on the Fourier series representation of /. Each approximation 
F/v,m is the quotient of a trigonometric polynomial of degree N and a trigonometric polynomial of degree 
M. It was proven that these “Fourier- Pade” approximations converge point-wise to (/(x + ) + f(x~))/ 2 more 
rapidly (in some cases by a factor of l/k 2M ) than the Fourier series partial sums on which they are based. 
Although these approximations do not “eliminate” Gibbs phenomena, they do mitigate its effect. In partic- 
ular, the asymptotic value of the magnitude of the overshoot is reduced to about 6%, and, outside a “small” 
neighborhood of a point of discontinuity of /, the “unwanted” oscillations can (for practical purposes) be 
eliminated. 

More recently, Geer and Banerjee [6] introduced a new, simple class of periodic “singular basis func- 
tions”, which have special “built-in” singularities. They prove that these functions can be used to construct 
a sequence of approximations which converges exponentially to / in the maximum norm. In particular, this 
implies that the Gibbs phenomena can be completely eliminated, even when / has several points of discon- 
tinuity in the interval [— 7r, tt]. In order to construct these approximations, the locations and magnitudes of 
the jumps in / and its derivatives must be known a priori. For many applications, this may present a major 
limitation on the practical implementation of their method. However, an extension of their analysis shows 
that, if sufficiently accurate approximations to the locations and magnitudes of the discontinuities of / can 
be obtained, then the singular basis functions can still be used to construct a sequence of approximations 
which converges exponentially to / in the maximum norm, outside the union of a finite number of small, 
open sub- intervals containing the points of singularity. 

The main purpose of this paper is to present and analyze a simple, accurate, and robust technique to 
estimate the points of singularity and the associated jumps in a Lipschitz function and its derivatives, given 
only a finite number of its Fourier coefficients. A number of other investigators, including Gottlieb, et. al. 
([7], [8], [9]), Solomonoff [21], Eckhoff ([3], [4]) and Bauer [1], have investigated this problem, using a variety 
of approaches. In later sections, we shall compare our methods and results with theirs. After presenting 
our technique, we will illustrate how our estimates can be used with the singular basis functions method to 
accurately reconstruct the original (discontinuous) function, as well as its first few derivatives. 

To fix notation, let / be a 27r-periodic Lipschitz function. Then / has enough smoothness and regularity 
properties so that its Fourier series exists and converges to (/(x + ) + f(x~))/ 2, for every x £ [— 7 r, 7 r], i.e., 

1 N 
toa, f n(x) = ^ (/( a;+ ) + f( x ~)) i F n (x) = y + a.j cos (jx) + bj sin (jx), 

j - 1 



Here /(x+) ( f(x~ )) denotes the limit of / from the right (left) at x, and Fn denotes the N th Fourier series 
partial sum associated with / (x) . 

We shall say that xo is a point of simple discontinuity of / if 

[/(lo)] = /(*o ) - f ( x o ) / °- 


Also, we shall say that xo is a point of contact discontinuity of order q of / if 

[f W ( x o)] = o, k = 0, 1, and [/ (<r) (x 0 )] * 0. 
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Here / ^ denotes the k th derivative of /, and, since we are assuming that / is Lipschitz, both the left and 
right limits of each derivative of / exist at every point in the interval [— 7 r, 7 r]. We call a point where / has 
either a simple discontinuity or a contact discontinuity a point of singularity of the function. 

We now assume that / has a finite (unknown) number, n, of singularities which lie at the (unknown) 
locations x\, X 2 , ... , x n in the interval (— 7r,7r]. The fundamental problem we shall address and solve in 
this paper is that of determining accurate and robust approximations of the singularity locations {x s } and 
the associated jumps [/^(x s )] , k = 0, 1 , given only the first 2 TV + 1 Fourier coefficients {dj,bj} . In 
Section 2, we show how initial estimates of these quantities can be obtained for the simple discontinuities 
of / (i.e., for k = 0), by capitalizing on certain properties of Gibbs phenomenon near these points. In 
Section 3, we formulate a least-squares optimization problem, which fits an asymptotic form of the Fourier 
coefficients, involving the points {x 5 } and {[/(x a )]} , to some of the known coefficients. It is shown that this 
optimization problem can be solved efficiently, using the results of Section 2 as initial estimates. The result 
of this least-squares fit provides estimates of the singularity parameters of significantly higher quality than 
those obtained without the least-squares approach. In Sections 4 and 5, these results are extended to finding 
the jumps in f'(x) at each point of simple discontinuity of /, as well as the locations and magnitudes of the 
first order contact discontinuities of /. In Section 6, these results are generalized to find discontinuities in 
the M th derivative of /, for M >2. We shall refer to our method, including the determination of both initial 
and improved estimates of the singularity parameters, as the “least- squares parameter estimation” (LSPE) 
method. We compare our results with those of other investigators in Section 7. In Section 8, some aspects 
of the robustness of the LSPE method are discussed and illustrated with three examples. In Section 9, it is 
shown how the LSPE method can be used with the singular basis functions introduced in [6] to reconstruct 
the original (discontinuous) function, as well as its first few derivatives, with exponential accuracy. Some of 
the details of our analysis are included in two appendices. 


2. Points of simple discontinuity - initial estimates. It is well known that, around a point, x s , of 
simple discontinuity of /, the Fourier series partial sums { F n] exhibit Gibbs phenomenon. This phenomenon 
includes the oscillatory behavior of Fj v near x s , as well as the tendency of Fn to “overshoot” the magnitude 
of the jump in / at x 3 . The difference between the extrema of F n closest to x s asymptotically approaches 
about 118% of [f(x s )] 1 as N — * oo. Also, the locations of the extrema of Fn (x) , nearest to x Sy occur at 
x — T] s ± = x s ± ir/(N + 1) + 0(iV -2 ), as N — ► oo. More precisely (see, e.g., Carslaw [2]) 

lim F n (rj 3<± ) = ifo +) + /<*«-) ± 1 Si(Tr) [/ (*,)] , 

TV — ►oo Z IX 

and hence 


( 2 . 1 ) 

Here 


lim {F n ( v a ,+ ) ~ Fn (v,-)} = r Si 00 1/ (*.)] • 


N—*oo 


2 000 = 2 f52» 

7T 7T J Q U 


7X 


du = 1.17898. 


( 2 . 2 ) 


We now observe that, if we define a new sequence of partial sums 

F n (x + n/(N + 1)) - F N (x — 7x/(N + 1)) 


Dn (^) — 


(2/tt) Si(7T) 


then, as N — > oo, Dn (x) — > 0, if / is continuous at x, but Dn (x) — > [f (x*)] , if x = x s is a point of simple 
discontinuity of /. Thus, a simple method of identifying initial estimates of the points of simple discontinuity 
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of /, as well as the corresponding jumps, is to find those points in the interval ( — 7r, 7 r] for which the sequence 
{Dn (a:)} converges to a non-zero constant, as N — ► oo. At these points, D n exhibits peaks, which become 
narrower, with increasing N. More specifically, initial estimates of the points of simple discontinuity to / 
can be obtained by locating the extrema, say x a , of .D/v, such that the corresponding values D^(x s ) do not 
tend to zero, as TV — ► oo. For such a point, D^(x a ) provides an estimate of [/ (x a )] . 

We now illustrate this idea with two examples. 


Example 1 . Let / be defined by 


(2.3) 


/(*) 


1 -x, 

< 5x 2 - 37x + 67, 

x 3 — 15x 2 + 75x — 125, 


0 < x < 1, 

1 < x < 3, 

3 < x < 4, 

4 < x < 5, 

5 < x < 27r, 


f(x + 2tt) = f(x). 


This function has only one point of simple discontinuity, namely, at x\ = 3. However, it has multiple contact 
discontinuities; one of order 1 at x 2 = 1, one of order 2 at x 3 — 4, and one of order 3 at x 4 = 5. The jumps 
in / and its derivatives at these points are summarized in Table 1. 


Points—* 

Jumpsi 

£2 — 1 

x\ = 3 

x 3 = 4 

£4 = 5 

[/(•)] 

0 

3 

0 

0 

[/'(■)] 

-1 

-6 

0 

0 

(/"(•)] 

0 

10 

-16 

0 

[/"'(•)] 

0 

0 

6 

-6 


Table Is Exact locations and magnitudes of the discontinuities for Example 1. 


A representative plot of D^(x) for this example is shown in Fig.L By finding the extremum of D N 
closest to x = 3, we can estimate both the location and magnitude of the simple discontinuity of / at this 
point. The absolute error in estimating the point of discontinuity, and the relative error in estimating the 
corresponding jump in /, using Dn, are summarized in Table 2 for several values of N. Also, in Figs. 2 and 
3, we have plotted these errors (suitably normalized) as a function of 1/N. These plots suggest that, using 
D n , we can approximate the location of the simple discontinuity of / with an error that is O (N " 2 ) , and 
that we can approximate the associated jump with an error that isO(A r_1 ),as7V— + oo. 


N 

Absolute Error in x\ 

Relative Error in D^(x 1 ) 

16 

3.2- icr 2 

3.48 • 10 -1 

32 

9.4- 10~ 3 

1.95- 10- 1 

64 

2.6- 10~ 3 

1.07 • 10- 1 

128 

6.8- icr 4 

5.34 • 10 -2 

256 

i.7- icr 4 

2.72 • 10~ 2 


Table 2: Errors in the initial approximations to the location and magnitude of the simple discontinuity in 

Example 1 using the partial sums {Dtv(x)} . 

Example 2. Let / be defined by 
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(2.4) 


/(*)=< 


0, 0 < x < 1, 

e x , 1 < x < 2, 

sin (|) , 2 < x < 5, 

0, 5 < x < 27 r, 


/(x + 2 tt) = /(x). 


Here, / has three points of singularity (with non-zero values of [/^(x)] at each of these points), and is 
analytic between these points. (This example has been used as a “benchmark” example by several other 
investigators. We will use it in later sections to compare our results with theirs.) The errors in locating 
the point of simple discontinuity at x\ = 1, along with the errors in approximating the corresponding jump 
[/(x i)] , are summarized in Table 3 for several values of N. The magnitude of the errors corresponding to 
the simple discontinuities at ~ 2 and X3 = 5 are similar to those reported in Table 3. When these errors 
are plotted as functions of 1/AT , we obtain graphs that are qualitatively very similar to those shown in Figs. 
2 and 3. 


It can be shown that the rates of convergence of our approximations observed in these two examples hold 
in a more general setting. In Appendix A, we outline a proof of the following theorem, which establishes 
rigorously the rates of convergence of our approximations to {x s } and {[/(x s )]}, as N — * 00, using the partial 
sums Dn(x). 


N 

Absolute Error in x\ 

Relative Error in Dn(x 1) 

16 

1.9 • icr 2 

1.8 • 10" 1 

32 

6.6- icr 3 

8.8 ■ 1(T 2 

64 

i.5- icr 3 

4.2 • 10“ 2 

128 

3.6 • 10- 4 

2.1 • 10“ 2 

256 

8.8 • 10" 5 

1.0- 10 -2 


Table 3: Errors in the initial approximations to the location and magnitude of the simple discontinuity at 

x\ = 1 in Example 2 using the partial sums {Z)^(x)} . 


THEOREM 2.1. Let f(x) be a 2n -periodic Lipschitz function with no points {x a } of simple discontinuity 
in the interval (— 7r, tt]. Then , as N — » 00, the partial sums Dn (x) , defined in Eq.(2.2), converge to zero 
for all x ^ x s , and converge to [f(x s )] at x = x s . The point x 8 , at which Dn (x) attains its extremum in a 
neighborhood of x 5 , approximates the point x s to within an error which is O ( N ~ 2 ) , and the corresponding 
value Dn(x 3 ) approximates [/(x a )] to within an error which is O ( N ~ l ) , as N — ► 00. 

3. Points of simple discontinuity - improved estimates. We assume that the number no of simple 
discontinuities of /, as well as initial estimates x s and Jo )S = Dn(x 8 ) of x s and [/(x 5 )] , respectively, for 
s = 1 , 2, ..., no, have been obtained using the method of the previous section. We now show how to improve 
the quality of these estimates, by utilizing the asymptotic form of the coefficients {a^, bj} , expressed in terms 
of the parameters {x s } and {[/(x s )]} , and a least-squares approximation technique. 

The proof of the following lemma follows from the definitions (1.1) and by the repeated use of integration 
by parts (see, e.g., Kreyszig [14], pp. 489-493). 

Lemma 3.1. Let f be a 2n-periodic Lipschitz function, with singularities at a finite number of points, say 
at x — x 3 , s — 1, ..., n, in the interval (—7 r, 7r]. Then, for any non-negative integer M , the Fourier coefficients 
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of f can be expressed as 


= \ E | sin (^*) ^ E { jlLi [/ (2fc) (^.)] j 


'L(M-1)/2J fc+1 

+ cos(jx s )[ Y1 ~j 2 kW [f {2k+1) M 


k = 0 


+ o(i/j M+2 ), 


(3.1) 


1 » ( /W2J ,_ lV t 

b i = ~ E \ cos (j x ‘) E TalfeTT [/ (2t> (*.) 

J >=1 1 \ k =0 J 


( L (M-1)/2J fc+1 _\1 


as j — * oo. In the upper limits of the inner summations, [tfj denotes the greatest integer not exceeding q. 
We use Lemma 3.1 with M — 0 to write 


jncij + ^ sin (jx s ) if (*«)] = ° (j ') . 


5=1 


n 0 

(3.2) jnbj - ^Pcos (jx„) [f (x,)] = O (j -1 ) , 

5=1 

as j — * oo, and then define a weighted least-squares error function 

N 

(3.3) E = E (xi,...,x no , Jo, no ) = E w i E j 0) > 

j = + 

C np \ 2 / no 

jiraj + s'm(jx s ) J 0 , s I + I jirb, - Y^cos (jx,)J 0ta 

Now, given the coefficients {a Jy bj}, we seek to determine the values of i\ 7 ... , x no , Jo,i> ••• , Jo,n 0 > sa y x \i 
... , x* 0 , J 0 * l5 ... , no , suc ^ that E is a minimum , i.e., 

(3.5) X no j */o,l> ***» *^0,n o l = (^li •••? ^no ? *^0,li ***> *^0,no) * 

& 1 >^0,T»Q 

We note that, in order to have a true least-squares minimization problem, we must require that R > no, 
since there are 2no parameters to be estimated. However, we must also require that N — R 1, in order for 
the asymptotic form of the Fourier coefficients to be valid. (We find the weighted least-squares strategy to 
be particularly useful when we assign larger weights to higher order terms defined in Eq.(3.3). This is due 
to the increasing accuracy of the asymptotic form (3.2) with the index j . ) 

From Eqs.(3.2)-(3.5), it follows that x* — ► x s and Jq 8 — + [f {x B )]> for s = l,...,no, as N — ► oo. More 
precisely, for the case uj = 1, it follows from these equations that 
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(3.6) 


oo. 


i: = I ‘ + FwT ^ +0<1/w4) ' “ d 

Jo,s = if (*.)] - + O (1 /N 3 ) , as N-* 

Thus, the parameters we wish to estimate have now been characterized as the solution to a certain 
optimization problem. We shall demonstrate in later chapters that this formulation offers some nice advan- 
tages over the approaches used by several other investigators. In particular, we can see, at least intuitively, 
that the approximations determined by this approach should be less sensitive to “small errors” in the coef- 
ficients {a-j^bj } , due both to the least-squares nature of the problem, and to the fact that the definition of 
E “averages” over several (i.e., 2 R) coefficients. 

To perform the minimization required by Eq.(3.5) in an efficient manner, we use the estimates ob- 
tained by the method of Section 2 as initial estimates for the nonlinear weighted least-squares method of 
Levenberg-Marquardt ([18], [19], [20]). Although the function in Eq.(3.5) can be minimized using a general 
unconstrained optimization technique, the special structure of the gradient, the Hessian matrix, and the gen- 
eral least-squares nature of E is exploited by the Levenberg-Marquardt method. In choosing an optimization 
method, it is important to note that E is a badly scaled objective function, which exhibits highly oscillatory 
behavior in some directions and very smooth, slowly changing behavior in other directions. In particular, a 
“flat” region of E creates a problem for derivative based methods, because of the closeness of the derivatives 
to zero in that region. To illustrate this behavior, in Fig. 4 we have plotted E for Example 1 as a function of 
x\ and J 0 ,i, near the minimum of E. In the figure we have set N — 32 and R — 12. For smaller values of R , 
the oscillations in the x\ direction are even more pronounced. The problem of the ill-conditioned associated 
Hessian matrix for E can be partially mitigated by assigning larger values to the weights Wj, with increasing 
index j, in Eq.(3.3). 

From these observations, it is important to note that we require good initial estimates of {x s } and 
{[/(x s )]} to initialize the search technique. Fortunately, the initial estimates obtained using the method 
based on the partial sums Dj y, described in Section 2, are good and, in fact, are even better than we actually 
require. For example, from the definition of E (see, also, Fig. 4) we see that the width of each “valley” in the 
ir s -direction is 0(1 /N). Thus, to begin a search in the “proper valley”, the initial estimate of each x s should 
lie within a distance that is at least 0(1/ N) of x s . From Theorem 2.1, we see that our initial estimates lie 
within a distance which is 0(1 /N 2 ) of x s . 

To illustrate this method, we again consider Examples 1 and 2. For Example 1, we use the initial 
estimates obtained in Section 2, whose errors are shown in Table 2, as starting values for the weighted 
least-squares method. The errors in the approximations to the discontinuity location x\ — 3, and the 
corresponding jump in /, are summarized in Table 4 for several values of N , and are plotted (with suitable 
normalization) in Figs. 2 and 3. In Eq.(3.3), each Wj = j. 
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N(R) 

Absolute Error in x* 

Relative Error in Jq x 

16 ( 8 ) 

1.14- 10" 3 

1.02 • 10~ 2 

32 (12) 

2.7- 10 “ 3 

2.45 • 10 - 3 

64 (15) 

6.05 • 10 ~ 4 

2.91 • 10 - 4 

128 ( 20 ) 

1.45- 10 - 4 

4.86 • 10 ~ 5 

256 (28) 

3.38 • 10 “ 5 

2.45 • 10 - 5 


Table 4: Errors in the estimates using the LSPE method of the location and magnitude of the simple 

discontinuity for Example 1 . 

Tables 4 and 2 , as well as Figs. 2 and 3, illustrate that our approximations to both the location of 
the discontinuity and the associated jump have improved. In particular, they illustrate that, using the 
optimization idea, the order of the approximation to x\ remains 0(N~ 2 ) ) but with a smaller proportionality 
constant. (Using Eqs.(ll) and Eq.(A.9) of Appendix A, we see that this constant is smaller by a factor of 
( 7 rSi( 7 r )) -1 = 0.172.) The order of the error in the associated jump in / has improved from O (A -1 ) to 
O ( N ~ 2 ) . 

For Example 2 , we again use the estimates from Section 2 , whose errors are given in Table 3 , to start the 
weighted least-squares method. The errors in the discontinuity location x\ = 1, and the corresponding jump, 
for different fixed values of TV, are summarized in Table 5. The order of the errors for the discontinuities at 
X 2 = 2 and X 3 = 5 are similar. 


N(R) 

Absolute Error in x\ 

Relative Error in Jq x 

32 (12) 

1.59 ■ 10 ~ 3 

2.93 • 10 " 3 

64 (15) 

3.1 • 10 “ 4 

5.96 • 10 -4 

128 ( 20 ) 

7.22 • 10 “ 5 

9.52 • 10 - 5 

256 (28) 

1.63 ■ 10 " 5 

3.04 • 10 “ 5 


Table 5: Errors in the estimates using the LSPE method of the location and magnitude of the simple 

discontinuity at X\ = 1 for Example 2. 


4. Discontinuities in /'(*) - initial estimates. To find initial estimates of the discontinuities in /', 

we are tempted to apply the method of Section 2 directly to the differentiated Fourier series partial sums. 
However, due to the presence of simple discontinuities in /, the differentiated Fourier series does not converge 
in a neighborhood of these points. Lanczos’ ([16], pp. 61 - 74) “sigma factor smoothing” provides one way 
to mitigate this problem, but not to overcome it. The modification of the Fourier coefficients by the sigma 
factors causes the differentiated series to converge at all points where /' (x) is finite. However, even after 
smoothing, the series diverges at all points where the derivative of the original function does not exist. 

We now show that a slight modification and extension of the method based on the partial sums {D^} , 
described in Section 2 , can give us initial estimates of the discontinuities in /' which are sufficiently accurate 
for our purposes. To see how this can be done, we use the (improved) estimates {x*, Jq a } obtained in the 
previous section to define a new sequence of Fourier coefficients j by 


(4.1) 


1 "o . no 

a > — a j + ^ Jo, s sin(j/x a ), b [ ^ = b j — — Jq s cos (jx*), 

— 1 3 S=1 


8 = 1 


for j = 1 , 2, ..., TV. (Intuitively, all we are doing here is “subtracting off” our best estimate of the parts of the 


8 












coefficients related to the simple discontinuities of /.) We then define a new sequence of partial sums 


d N 

F\ ,N = ^ |aj 0) cos(jx) + 6j 0) sin(jx) j 


j=i 


(4.2) 


N 

= Y cos(jx) + b { p sin(jx)} , 


3 = 1 


(1) _ -i(0) , (1) _ . (0) 

a ) \ b ) = -J a j 1 

which we can think of as the derivative of our original partial sum, but with the effect of the simple discon- 
tinuities “subtracted off 1 ’ . The coefficients in have the asymptotic form 


1 ( _ n ° n o \ 

a i 1} = - ( X] cos U x ») if ( x »)] - 5Z c ° s U x s) j o,s ) 

\ S~ 1 8 = 1 ) 


(4.3) 


sin (J x * ) [/' ( x *)l + ° (j 2 ) - 

J s=l 


+ / n o no 

h( P = - f sin (i x s)lf (**)] - X^sin (jx*)J£ s 


\s=l 


S— 1 


+ ~^ 2 C0S ^ X ^ [/'( I *)l + °(i 2 )> 

*' 5 = 1 

as j — > oo. Here n\ is the (unknown) number of points of discontinuity in /'. (We shall assume that n\ > no, 
by including all of the points of simple discontinuity of / in the set of points where /' may be discontinuous.) 
Using Eqs.(ll), we write Eqs.(4.3) as 


a T = l 5/' ( X »)l (^2 - j) Sin (j' x ») + ° (772) 

1 "I 

i/' (^s)] sin (j^ s ) +o(j~ 2 ), 

J »=no + l 

( 4-4 ) h T = l (*•)) (j - cos u x s) + o (^) 

+= (/' ( x *)] cos(jx a ) +0(j - 2 ), 

5 — n 0 f 1 


as j — > 00 . (Here the points {x 3 } , s = no + 1, ni, correspond to contact singularities of / of order one.) 

We now make two observations. First, by comparing Eqs.(4.4) and Eqs.(3.1), we see that, for x near 
the contact singularity at x s , s = ra 0 + 1, ...,n 1 , the partial sums {Fi yN } behave like a function which has a 
simple discontinuity at this point of magnitude [/' (x a )]. Thus, if we define a new partial sum Di y ^(x) as in 
Eq.(2.2), but based on the partial sum F\ y N, instead of F/v, i.e., 


(4.5) 


D lyN ( x ) 


F\,n {x + k/{N + 1 )) - F lyN (x - 7 r/(7V + 1)) 
(2/tt) Si( 7 r) 


9 



then it follows that 


Di,n{x s ) — > [/'(a:,,)] , at a point of contact discontinuity, 


as N — ► oo. Second, to investigate the behavior of D\^(x) near the point of simple discontinuity at x s , 
s = 1, 2, ..., no, we first observe (following the same reasoning as in Carslaw [2]) that 


N 


lim F un (t) s ,±) = 


/'(*+)+/'(* 7) 


7T \ 7 T / 


Thus, it follows that 




at a point of simple discontinuity, 


as N — ► oo. (Here (Si(7r) — l/7r)/ Si(7r) = 0.82812.) For any point x that is not a point of simple discontinuity 
or contact discontinuity of /, Di^n(x) — * 0, as TV — ► oo. 

To illustrate these ideas, we use them to find initial estimates of points of first order singularity of / for 
Example 1. In Fig. 5, we have plotted D i ^ for this example with several values of N. The absolute errors in 
the estimate of the point of contact discontinuity at = 1 , and the relative errors in the estimation of the 

associated jump [/'(l)] , using the partial sums of are summarized in Table 6 for several values of N. 

The relative errors in the estimation of the jump [/' (3)] are also listed. Similar to the method of Section 
2, we can show, in general, that, using the partial sums {D \^} , the error in approximating the locations 
of the points of contact discontinuities is 0(N~ 2 ), while the error in the approximations to all of the jumps 
[/'(*.)] ^ 0(N~'). 


N 

Abs. Error in X2 

Rel. Error 

in -Di,/v(x 2 )] 

Rel. Error in f?i,jv(xi)] 

16 

2.1 • 10~ 3 

7.4 

1(T 3 

4.5- 10" 1 

32 

2.4 lO -5 

1.5 

10“ 3 

2.2 ■ lO" 1 

64 

1.2 10- 5 

7.8 

• 1(T 5 

1.1 ■ lO" 1 

128 

5.4 lO -7 

2.0 

• 10“ 6 

5.1 • 10“ 2 

256 

2.7- 10~ 8 

2.9 

• io - 6 

2.5 • 10 -2 


Table 6: Errors in the initial estimates of the discontinuities in /' for Example 1 using the partial sums 

{Di,n} • 


5. Discontinuities in f'(x) - improved estimates. Once the number of discontinuities in / and /' 
are determined, as well as initial estimates of the locations and magnitudes of all of the jumps in / and /', we 
again seek to improve these estimates by using a least-squares fit of the asymptotic form of the coefficients 
to the given coefficients. Using Lemma 3.1 with M = 1, we now define 


(5.1) 


where 


N 

~ E ^Xi , X ni , t/o,!, Jo,ni j «A,lj *A,ni ^ ^ \ 


j=N+l-R 


= ^j 2 naj + (i sm(jx s )J 0 , s + cos(jx s ) + 


10 




(5.2) 


2 


' f 

j nbj - 22 jjcos (jx s )J 0 ,s - sin (jx a )J hs 


8=1 


(Note that, for ease of notation, we have let the index s range up to n i, with the understanding that = 0, 
for s = no + 1, ni.) Then, given the coefficients {flj, bj} , we seek to determine the values of jx s , Jo, s , Ji ?5 1 , 
say {x*, Jq s) J* 5 } , s = 1, 2, such that E is a minimum , i.e., 


(5.3) 


min E \ 


aji , Jo.i , j *7] 


, ni ) J* ni ). 


(Since there are now 3m parameters to be estimated, we require R > 3ni/2.) 
From Eqs.(5.1)-(5.3), it follows, for the case when uij = 1, that 


x* s — x s + O (l /N 5 ) , 1 < s < n 0 , 


x* s = x s + O ( l/N 3 ) , no + 1 < s < m ? 

(5-4) ^ 0 , s = [/(^)]-[/"(a: s )]^+O(l/iV 3 ), 1 < « < n 0 , 

Jl a = [f'(x s )} + 0(l/N 2 ), l<s<n u 


as N —> oo. 

In Table 7, the absolute error in estimating the point of contact discontinuity at X 2 = 1 
and the relative errors in the estimation of the jumps [/'( 1)] and [/'(3)] , using the method 
are summarized. 


N(R) 

x\ 

T* 

J l ,2 


J* 

J 0,l 

Ki 

32 (15) 

2.1 • 10~ 4 

1.47 - 10“ 2 

3.29 • 10“ 5 

4.16 - 10“ 3 

1.97 - 10“ 2 

64 (18) 

8.14 • 10~ 5 

2.0 -10- 3 

1.95- 10“ 6 

CO 

1 

O 

r-H 

0 

q 

1— 1 

2.72 • 10“ 3 

128 (20) 

1.93- 10- 6 

7.19- 10“ 4 

7.24 • 10- 7 

to 

t— * 

0 

l 

3.17 - 10- 3 

256 (46) 

6.4- 10~ 7 

1.56- 10“ 4 

2.00 • 10“ 9 

6.41 • 10- 5 

2.17 • 10“ 4 


Table 7: Errors in the estimates of the discontinuities in / and /' for Example 1 using the LSPE method 

with M = 1. 

Comparing Table 7 to Table 6 illustrates that the errors in the estimates of [f f (x 1 )] are significantly 
smaller using the least-squares approach. Of course, we also obtain “new” estimates of the location and 
magnitude of the discontinuity in / at x\ = 3. The corresponding errors are also shown in Table 7. We 
observe that the errors in the estimates of [f f (x 2 )], and the location of the first order contact discontinuity 
(at X 2 = l) are somewhat greater when compared with the errors in the initial estimates they use. This 
should not, however, be a cause of concern, as the estimates that are obtained using the LSPE method are 
good y and are within the order of accuracy that is predicted for them (see Eqs.(5.4)). 

6. Estimation of discontinuities of order M, for M > 2. Assume now that we have used the 
methods of the previous sections and that we have obtained, by the LSPE method, estimates {x*} and 
the locations {x s } and the associated jumps [/^(x s )] , respectively, of the singularities of /, for 


for Example 1, 
just described, 


ll 




k = 0,1, ..., M l, and 5 = 1,2, We define 


■ L 


'L(M-1)/2| . 


k = 0 


([(M-V/ll k+1 \ 

5 Ux») ( Z j2fc+2 ^fc+l.a 1 | 


and the partial sums 


where 


n M -i /L(M-l)/2j fc 

& r j) = 6 ; - - z | c ° s (i x s) ( z pkTT j k. 

S=1 V fc=0 7 


+ sin(ja:* ) Z 


fc=0 

'[(A/-2)/2j , ,, fc+1 


(_l)fe+l 
,-2fc+2 ,> 2fc+l,« 


k—O 


/ I \ iM iV 

Fm,n(x) = Z { a j M X) cospx) + b\ M ~ l> sin(jx)| 

^ ' j=i 

N 

= Z { a J M) cos 0' x ) + 6 j M) sin (j x )} . 


a (M) = ( _ 1)L M/2j 6 (M) ^ ( _ 1) LM/2J i M 6 (M-D | M eveQ) 

a (M) ^ ( _ 1)L M/2J 6 (M) ^ m odd 

Here, in general, rij denotes the number of points in the interval ( — tt, 7r] where has a discontinuity for 
at least one value of A: < j. 

To find initial estimates of the locations and magnitudes of the discontinuities in we define 

n x Fm,n(x + k/(N+1))- F m ,n{x-tt/(N +1)) 

Dmn W ” (2/»)SiM • 

Then it follows that 


Dm,n(Xs) 


(*m,s ■ [/ (M) ( X »)] , S = l,2,...,n M , 
0, otherwise, 


as N — ► oo. Here cum, a is a scalar that depends on both M and s . In particular, 


and 


&2 f s 


(Si(7r) - 1 /n)/ Si(7r), 1 < s < no, 

1, no + 1 < s < n 2 , 


^3, a 


' (Si(7r) — l/7r — 6/ 7r 3 )/ Si(?r), 
< (Si(7r) - 2/n)/ Si(ir), 

1, 


1 < s < no, 
no + 1 < s < ni, 

7l\ + 1 < S < 71 3. 
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Thus, Dm,n{%) can be used to find initial estimates of the locations and magnitudes of the jumps in /( M )(:r), 
as in Sections 2 and 4. 

To illustrate these ideas, we set M = 2 and use D 2 ,n(x) to estimate the locations and magnitudes of 
the jumps in /" for Example 1. The results are summarized in Table 8. (We find that D 2 % n( x ) — ► 0 in a 
neighborhood of x = x 2 = 1, suggesting that [}"{x 2 )] is zero.) 


N 

Abs. Error 

in D 2 ,n{x 2 ) 

Abs. Error in X3 

Rel. Error in D 2 ,n{x 3) 

32 

2.57 

lfT 1 

1.94 • 1(T 3 

3.07 x 10“ 2 

64 

1.03 

10- 1 

5.07 x 10~ 4 

1.53 x 10 -2 

128 

3.49 

• 10- 1 

1.33 x 10“ 4 

7.73 x 10“ 3 

256 

3.03 

10- 1 

3.30 • 10“ 5 

3.89 • 10~ 3 


Table 8: Errors in the initial estimate of the discontinuities in /" for Example 1 using the partial sums 

{£>2, at} • 


Once initial estimates of the locations and magnitudes of the jumps in have been determined, we 
again improve these estimates by using a least-squares fit to the asymptotic form of the coefficients. Thus, 
we define 


( 6 . 1 ) 


N 

^ E 

j=N-R+l 


“jEj 


(M) 


( 6 . 2 ) 


E. 


(M) 


n\f 

j M+l *aj - J2 

5 = 1 


{ [M/2] 

E (~ 1 ) k+1 3 M ~ 2k J2k,s } + 
k=o 

-1 \ 2 


cos 


[(M— l)/2] 

(jx>) { E (- 1 ) fc+1 i M ' 1_2fc ^2fc+l, i 


0 

tim 


+ 


•M + l 


nb 3 - E 


5=1 


[M/2] 


COS (jx s ) < ^2 2k J2k,s > + 


k = 0 


-| \ 2 


{ [(Af — 1)2] 

E (- l ) fc+1 j M_1_2fc i 2fc+li , 

k = 0 

and choose values of the n M { 2 + M) parameters {x s } and j , say, |x*, , s - 1,2 k 

0, 1, so that E is a minimum , i.e., 

min E (xi , Jo.i, JM,n M ^ = E (aq , Jq i, • 

*4, a X 7 


Here i? > hm(2 + M)/ 2, with iV - i? » 1. 

When we apply these ideas to Example 1 with M = 2, we obtain the results whose errors are summarized 
in Tables 9 and 10. (We find that the estimates obtained for [/"(a^)] suggest, rightfully, that [/"(x 2 )] = 0.) 

Setting M — 3, we find initial and improved estimates of the contact singularity of order 3 at X 4 = 5 
for Example 1. In this case, the maximum error for the location or associated jump in any of the improved 
estimates is of the order of 10“ 13 . To understand why the errors should be so small, we note that, since / 


13 




of Example 1 is piece-wise polynomial, of degree at most 3, the assumed form of the Fourier coefficients in 
Eq.(6.2) is exact for M > 3. Hence, the estimates obtained by the method are “essentially exact”, to within 
the error inherent in the minimization procedure. 


■ 

x\ 

7* 

J 0,l 

J 1.1 

T* 

J 2,l 



2.92 • 10- 5 

5.33 • 10- 4 

5.4 • 10“ 3 



1.33 • 10~ 6 

1.03 • 10~ 4 

3.22 10- 3 

128(32) 

1.0- 1(T 9 

1.03 • 10“ 7 

5.67 • 10" 6 

9.8 • 10“ 4 

256(56) 

8.75- 10- 11 

2.73- 10“ 8 

1.67 ■ 10~ 6 

4.5 ■ 10“ 4 


Tkble 9: Errors in the estimates of the discontinuity at x\ =3 and the associated jumps for Example 1 

using the LSPE method with M = 2. 


N(R) 

x ' 2 

7* 

J l,2 

r* 

*“'2,2 

*3 

7* 

J 2,3 

32(18) 

3.95 • 10“ 5 

5.72 - 10" 4 

2.3 • 10~ 2 

5.4 • 10“ 4 

7.75 • io~ 4 

64(22) 

3.0 • 10“ 6 

8.09 • 10- 6 

9.1 • 10~ 3 

1.14- 10~ 4 

1.99 • 10- 4 

128(32) 

5.0- IQ" 7 

2.8 ■ lO" 6 

6.67 ■ 10~ 3 

2.9 - lO" 5 

1.25 • 10“ 4 



1.5 • lO" 6 

1.43 • lO -3 

7.06 • 10 6 

1.24- 10" 5 


Table 10: Errors in the estimates of the discontinuities at X 2 = 1 and X 3 = 4, and the associated jumps, 

for Example 1 using the LSPE method with M = 2. 


N{R) 

x\ 

1 * 

JU 

32(15) 

3.29 - 10“ 5 

1.55 • 10- 3 

2.44 • 10 -2 

64(18) 

1.17- 10~ 6 

3.49 • 10- 4 

3.28 • 10“ 3 

128(24) 

1.65 - 10“ 8 

7.35 • 10- 5 

2.65 • 10“ 4 

256(41) 

2.15- 10“ 9 

1.77 - 10“ 5 

1.63 • 10“ 4 


Table 11: Errors in the estimates of the discontinuity at X\ — 1 for Example 2 using the LSPE method 

with M — 1. 


On applying the above mentioned ideas to Example 2 with M = 1 and M = 2, we obtain results that 
are qualitatively similar to those for Example 1. The results for the discontinuity at X\ = 1 are summarized 
in Tables 11 and 12. The results for the other points of discontinuity in / are qualitatively similar to those 
reported in Tables 11 and 12. 


N(R) 

x\ 

7* 

J 0,l 


7 * 

^2,1 

32(15) 

2.96 ■ 10“ 6 

1.29 - 10“ 5 

3.54 • 10“ 3 

1.38 • 10“ 2 

64(22) 

1.28- 10~ 7 

8.29 - 10“ 7 

7.29 • lO" 4 

3.34 • 10“ 3 


6.47 - 10" 9 

2.26 • lO" 8 

1.63 • 10" 4 

3.94 • lO" 4 


3.72 • 10“ 10 

2.92 - 10“ 9 

3.91 • 10“ 5 

1.16 ■ 10~ 4 


Table 12: Errors in the estimates of the discontinuity at X\ = 1 for Example 2 using the LSPE method 

with M = 2. 


For a general value of M > 0, we conjecture that the estimates obtained from the LSPE method satisfy 

x’ = x, + 0(N~ m ~ 2 ), a = l,2, —,n, 


n,. = [/ (fc) (x a )] +o(N- M ~ i+k ), 


(6.3) 
and 

(6.4) 
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k = 0,1,..., Af, 5 = 1,2, ...,n. 







































Although we do not have an analytical proof of these results, they have been verified using Mathematica 
[22] for 0 < M < 4 and for 1 < n < 6. Also, the results of several numerical experiments are consistent with 
them. 


7. Comparison with other methods. Three recent methods with a similar goal of obtaining esti- 
mates of the locations of the discontinuities of / from a finite number of its Fourier coefficients have been 
proposed by Eckhoff [3], Kvernadze [15], and Bauer [1], 

Eckhoff [3] develops an algebraic equation of degree n for the n singularity locations in each period of /. 
The coefficients in his equation are obtained by solving an algebraic system of n equations, determined by 
the known coefficients in the truncated Fourier series. If discontinuities in the derivatives of / are considered, 
in addition to the simple discontinuities of /, that algebraic system is nonlinear in the unknown parameters. 
The degree of the algebraic system depends on the desired order M for the reconstruction, with a higher 
value of M normally leading to a more accurate determination of the singularity locations. The jumps in / 
and its derivatives, up to the order M, can be obtained by solving an additional linear algebraic system of 
equations. However, the robustness of the method deteriorates with increasing values of M. This is due to 
the ill-conditioned nature of the equations that must be solved. 

Kvernadze [15] proposes an algorithm to determine the discontinuities and the corresponding jumps 
in / using certain identities based on the partial sums of its differentiated Fourier series. Kvernadze first 
analyzes an expansion formula for the approximation of a 27r-periodic, piecewise smooth function with one 
discontinuity. An appropriate linear combination of certain identities, obtained via derivatives of different 
orders, is then used to significantly improve the accuracy of the estimation. It is then suggested that 
Richardson’s extrapolation method be used to refine the accuracy even further. For a function with multiple 
discontinuities, Kvernadze establishes a formula that eliminates all but one discontinuity in the function, 
and then treats the new function as one with a single point of discontinuity. Kvernadze does not address 
aspects of the numerical complexity and robustness of his algorithm. 

Bauer [1] uses the idea of band pass filters to find the discontinuity locations. He makes no effort to 
estimate the associated jumps or the points of contact discontinuities. In his work, Bauer introduces the idea 
of a global filter and a local sub-cell filter. There is a trade off between these two filters. Smaller errors can be 
obtained using a local sub-cell filter, but the accuracy decreases if there is a contact discontinuity. His global 
filter can be computed “once” and stored in memory, but the local filer cannot be computed until initial 
estimates of the locations of the discontinuities are determined. A comparison of Eckhoff ’s, Kvernadze’s, 
and Bauer’s methods show that, for a given value of AT, Eckhoff ’s results are consistently more accurate 
and robust. Also, Bauer’s local- filter method appears to be unable to control the error if the function has a 
contact discontinuity. 

The LSPE method we are proposing is based on a simple idea, is remarkably robust (see the next 
section), and requires only “standard” optimization techniques to implement. The method appears to 
provide estimates of the discontinuity locations which are of the same order of accuracy (or better) than 
those of Eckhoff, Kvernadze, and Bauer. In addition, the LSPE method provides estimates (of high accuracy) 
of the associated jumps in / and its derivatives at the points of discontinuity. 

To illustrate this comparison, the function / of our Example 2 was also considered by Eckhoff, Kvernadze, 
and Bauer. The absolute value of the largest error in the estimation of the points of simple discontinuity of 
/, obtained by Bauer, by Eckhoff, by Kvernadze, and by the LSPE method are summarized in Tables 13-16, 
for different values of N and M. 


15 



N 

32 

64 

128 

256 

Local Filter 

i.4- ltr 5 



3.14- icr 9 

Global Filter 


8.77 • 10~ 5 




Table 13: Maximum errors in the estimates of the singularity locations for Example 2 using Bauer’s 

method. 


TV 

32 

64 

128 

256 

M = 0 

1 . 4 - 10" 3 

i 

o 

fH 

1— 1 
CO 

7.3 • 10“ 5 

1 . 8 - 10~ 5 

M = 1 

2.2 • 10~ 5 

1 . 0 - 10~ 6 

5.6 • 10“ 8 

2.7 - 10- 9 

M = 2 

4.0 • 10~ 6 

1 . 2 - 10- 7 

5.1 • 10~ 9 

2.7 - 10- 10 

M = 3 

1.2 - 10- 7 

2.8 • 10~ 10 

5.3 • lO " 11 * 

4.4 • IQ ” 5 


Table 14: Maximum errors in the estimates of the singularity locations for Example 2 using Eckhoff’s 

method. 


Tables 13-16 illustrate the convergence of all the methods. The results obtained by EckhofF’s method 
and the LSPE method are generally more accurate (for M > 2) than Bauer’s results and Kvernadze’s results. 
The robustness of Eckhoff’s method deteriorates when M > 3. This pattern becomes evident on comparison 
of results of his method to our results when M = 3. We find that for larger values for TV, Eckhoff’s method 
deteriorates significantly, while the LSPE method appears to show no such adverse effect. 


TV 

32 

64 

128 

256 



6.1 ■ 10“ 7 

1.4 • 10~ 8 

3.5 • 10 -11 


Table 15: Maximum errors in the estimates of the singularity locations for Example 2 using Kvernadze’s 

method. 


TV 

32 

64 

128 

256 

f||f]||j|J 


3.31 • 10~ 4 

7.73 • IQ" 5 

1.78 - 10“ 5 

2339 





M = 2 

2.95 • 10~ 6 

1.38 - 10“ 7 

7.1 • 10- 9 

4.26 - 10“ 10 

EB3 


4.14- io~ 9 

4.9 • 10" 11 

1.47 - 10~ 13 


Table 16: Maximum errors in the estimates of the singularity locations for Example 2 using the LSPE 

method. 

8. Robustness of the LSPE method. Intuitively, the LSPE method should possess good robustness 

characteristics, due primarily to the underlying least-squares nature of the method and the fact that the 
objective function is an “average” involving several Fourier coefficients. Although we have not carried out a 
detailed analysis of the robustness of the LSPE method, in this section we consider some of its robustness 
properties suggested by three more examples. 

Example 3. We consider again the function of Example 2, but we now contaminate its Fourier co- 
efficients by introducing some random errors . We then apply the LSPE method using these contaminated 
coefficients, and compare the results with those obtained using the original (“exact”) coefficients. Letting 

{oj, bj} denote the exact Fourier coefficients of /, we define the new coefficients 

j = a j “b e a,j i bj = bj + Cb,j > 1 < j ^ TV, 

where e a j and are independent, uniformly distributed random variables on the interval [— e, e], with e > 0 
specified. 
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The relative errors in the estimates of x\ = 1, as well as the corresponding jump in /, using the LSPE 
method with M = 0, N — 64, and R = 32, are summarized in Table 17, for a range of values of e. The 
estimates of x\ and [/ (rr i )] appear to improve only marginally with decreasing values of e in the range 
considered. 


e 

x\ 

7 * 

J 0,l 

10“ 2 

1.48 • 10~ 3 

CN 

1 

o 

i-H 

q 

10" 3 

2.0- 1(T 4 

3.723 10“ 3 

i<r 4 

3.73- lO" 4 

5.27- 10“ 4 

0 

3.92 • 1CT 4 

1.77 • 10~ 4 


Table 17: Errors in the estimates of the location and the magnitude of the discontinuity at x\ = 1 for 
Example 3 using contaminated coefficients in the LSPE method with M = 0. 


Results for the case M = 1, N = 64, and i? = 32, are summarized in Table 18. The estimates of X\ and 
[/ (xi)] appear to be relatively insensitive to the increase in the order for M, for the range of values of e 
considered. However, there is a noticeable improvement in the estimates of [f'(x i)] with decreasing e. This 
sensitivity of the higher order jump estimates to e is explained by considering the form of Eq.(5.1) that is 
used for the least-squares parameter estimation. The effect of a multiplier that is a power of j is likely to 
have an adverse effect on parameters that are dependent on higher precision digits in the Fourier coefficients. 
As a result, we see a noticeable increase in the accuracy of the estimate of [f f (x i)] with decreasing values of 
e. (Errors for the other singularities of / are comparable to those in Tables 17 and 18.) 


6 

x\ 

4),1 

Jil 

10“ 2 

7.52 ■ 10“ 3 

1.03- 

10“ 2 

12.31 



1.53- 

10 -3 

0.97 

10 -4 

6.92 • 10~ 5 

2.05* 

10“ 4 

0.1 

10“ 5 

8.69 • 10- 6 

3.72* 

1(T 4 

1.49 ■ 10“ 2 

0 

1.97- icr 6 

3.97* 

10“ 4 

5.51 - 1(T 3 


Table 18: Errors in the estimates of the discontinuity location and the associated jumps for Example 3 
using contaminated coefficients in the LSPE method with M = 1 . 

Example 4. To help assess the effect of the closeness of two points of discontinuity on the computed 
results, for 0 < a < 2 tt — 1 we define 


( 8 . 1 ) 


f(x) = 


x, 1 < x < 1 + a 
0, elsewhere in [0, 27 t] 


/(x + 2tt) = f(x). 


For this example / has discontinuities at x\ = 1 and at x 2 = 1 + a. The relative errors in the estimates 
of xi and x 2 , as well as the corresponding jumps in /, using the LSPE method with M — 0, TV = 64 and 
R — 15, are summarized in Table 19, for a range of values of a. All of the estimates of the parameters shown 
appear to be relatively insensitive to the value of a in the range considered. Results obtained by Bauer’s 
method, on the other hand, appear to be significantly more sensitive to the distance between two consecutive 
singularities, especially as a — ► 0. 
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a 


7* 

J l,0 

*2 

7* 

J 2,0 


3.8 ■ 10" 4 

6.89- 10“ 4 

3.14 • 10" 4 

2.97 • 10 -4 


3.1 • 1(T 4 

6.46- 10- 4 

1.38 • 10 -4 

2.58 • 10" 4 

0.9 

3.09 • KT 4 

3.68- 10" 4 

8.57 • 10" 5 

6.94 • 10~ 5 



1.91 • 10 -4 

3.37 • 10" 5 

1.31 • 10~ 4 


Table 19: Errors in the estimates of the location and magnitude of the discontinuities for Example 4 for a 
range of values of a using the LSPE method with M = 0. 

Example 5. We now illustrate some of the robustness of the LSPE method with respect to noise 
in sampled data . In this case, the Fourier coefficients of a function / are computed using a Fast Fourier 
Transform (FFT) method applied to a data set of (slightly erroneous) functional values at evenly spaced 
values of the independent variable . 

We define / as the “usual” step function, but with small random errors, as 


( 8 . 2 ) 


/(*) 


— 1 -h e i , 0 < x < 7r, 

1 + e 2 , n < x < 27r, 


/(x + 2tt) = f(x). 


Here ei and e 2 are independent random variables that are uniformly distributed over the interval [— e, e] , with 
e > 0 specified. Then / has two simple discontinuities, at x\ = n and at x 2 = 2n. We observe that there are 
two major sources of error in the computed Fourier coefficients. One source of error is the noise amplitude 
e, which is introduced into / and, hence, into the sampled functional values. Another source of error is the 
use of the FFT, which uses only a finite number of sample points. For example, using 2 N + 1 sample points, 
only the first 2N + 1 Fourier coefficients can be estimated, before abasing occurs (see, e.g., Hamming [13]). 
When error is introduced into the data, it is natural to expect the relative error in the computed Fourier 
coefficients to increase with increasing index. As a result, instead of using the last R coefficients, we find 
that it is better to base the LSPE method on some of the lower order coefficients. For our illustration, we 
used a sample of size 65 (from which we can estimate the first 32 {aj} and {bj}), but we only use the ninth 
through the fifteenth coefficients in the definition of E. The relative errors in the estimates of x\ and [f(x i)] 
using the LSPE method, with M = 0, are summarized in Table 20. 


N = 15, R = 7 

Abs. Error in x* 

Rel. Error in Jq x 

Exact, e = 0 

1.0 ■ 10- 17 

1.0- 10- 17 

FFT, e = 0 

1.33 - 10- 2 

5.5 • 10“ 3 

FFT, e = 10~ 2 

1.28- 10- 2 

7.5 • IQ" 3 


Table 20: Errors in the estimates of the location and the magnitude of the discontinuity for Example 5 
using the LSPE method with M = 0. The Fourier coefficients are computed from complete and sampled 

data. 


From Table 20, we note that there is virtually no error if the exact Fourier coefficients (e = 0) are used 
in the definition of E. If we compute the Fourier coefficients using the FFT based on the “exact” (e = 0) 
sampled data, we find that the relative error in the estimates of the point of simple discontinuity and the 
associated jump at x\ are of the order of 10 -2 and 10 -3 , respectively. However, even when we set e = 0.01 
and again use the FFT to compute the coefficients, we find the errors remain essentially unchanged from the 
case when e = 0. 
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9. Reconstruction of /. In a companion paper [6], we introduced and studied a class of 27r-periodic, 
singular basis functions, which have special “built-in 1 ’ discontinuities. In [6] it was proven that these func- 
tions can be used to construct a sequence of approximations to a discontinuous function, which converges 
exponentially to / in the maximum norm. However, the construction procedure requires knowledge of the 
exact locations and magnitudes of all of the discontinuities in / and its derivatives. 

In this section we briefly summarize how the approximating sequence should be modified when only 
estimates of locations and magnitudes of the discontinuities in / and its derivatives are known. In Appendix 
B, we outline how the main proof in [6] can be modified to show that the estimates we have obtained by 
the LSPE method are of sufficiently high quality so that we again obtain an approximating sequence which 
converges exponentially to / in the maximum norm, for x in the domain D. Here D consists of the interval 
[— 7r, 7r], with the union of certain “small” open intervals surrounding the points of discontinuity of / removed. 
(The measure of £) converges exponentially to 27r; see Appendix B for details.) In addition, the derivatives 
of this sequence converge exponentially to the corresponding derivatives of /, for x € D. 

The singular basis functions {S n (x)} are defined by 

s 2 k(x) = sin(x) (1 - cos(a;)) A:_1/2 = ^ b 2kJ sin (jx),, 

°c 

( 91 ) S 2k +i(x) = 7 2 ^ + Y)j 0 ~ cos(x)) fc+1/2 = 2fc ^ 1 ’ - + XI a 2 k+i,j cos [jx), 

1 

4/c+l J 

a 2k+l f j = (~l)* +1 , b 2 h,j = — 1 

n(4j 2 -(2 z+l) 2 ) 

»= 0 

for j — 0,1,2,... , and k = 0,1,2,... . (For convenience in some of the formulas below, we also define 
a 2 k,j = ^2 fc-j-i = 0, for k > 0.) It is straightforward to show that S m (x) is C Tn ~ 1 [— 7r, 7r], while the jump in 
its m th derivative at x = 0 is 1. 

Now let M be a nonnegative integer and let the 2n - periodic function / have possible discontinuities in 
f^ k \ for 0 < k < M, at x = x 5 , s = 1, 2, ..., n, where —n < x 3 <1 r. We define 

M n 

(9-2) S- M (x) 

fc=0 s=l 

where the constants {A£ are determined recursively by 

k — l 

(9.3) 4k,. = [ 5 i fc) (°)]> s = l,2,...n, k = 

i=0 

Here x* and J£ s are the estimates of x s and [/^(x s )] , respectively, determined by the LSPE method, 
as in sections 3, 5, and 6. (We note that, if x* 3 — x s and if J^ s = [/^( x *)] 5 then f(x) - S^(x) will be 
C M [-Tr,Tr], at least, and hence its Fourier series will converge at a faster rate than the Fourier series of /. 
See [6] for details.) Once S^(x) has been determined, we define the family of approximations / ^ N to / by 

(M) N 

(9.4) ImM x ) = S*m( x ) + a J M) cos U x ) + b { j M) sin (jx), 
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(9.5) 


M n 

•r=*i-EE A k,s { a k,j cos (jx*) - bk,j sin(jx*)}, 
fe=0 »=1 


M n 

b \ M) = sin(jx*) + b k<j cos(jx*)} , 

fc=0 a = l 

j = 0, 1, 2, TV . 

To illustrate the reconstruction method described above, we reconsider Example 1. We reconstruct / 
and a few of its derivatives using the first 2 TV 4- 1 of its Fourier coefficients. The estimates of the singularity 
locations and the associated jump parameters used are obtained by the LSPE method, with M = 2 and 
TV = 32. Figure 6 illustrates the excellent agreement between / (solid fine) and /2 32 (dashed fine) . From 
Figs. 7 and 8, it is evident that the first two derivatives of /2,32 also agree, to within the plotting accuracy, 
with the corresponding derivatives of /. Figure 8 illustrates that there is a slight deterioration in the level 
of agreement between /" (solid line) and (dashed line) at points that are close to the points of 

singularity of /". However, this deterioration is expected and can easily be eliminated by increasing M 
and/or TV. 

10. Conclusions, discussion, and future directions. A simple, accurate, and robust method (the 
LSPE method) has been introduced and studied to estimate the locations and magnitudes of the jumps in a 
function / and its derivatives, using only the information contained in the first 2 TV + 1 Fourier coefficients of 
/. These estimates can then be used with a simple class of periodic “singular basis functions” to construct a 
sequence of approximations which converges exponentially to /, and its derivatives, as TV — > oo, in the maxi- 
mum norm, outside the union of a finite number of small open intervals that contain the points of singularity 
of /. The total measure of the union of all such “small” intervals approaches zero exponentially as TV — ► oo. 
In particular, this implies that the effects of Gibbs phenomena can be completely eliminated , even when / 
has several points of discontinuity. Also, the singularities of / may be either points of simple discontinuity, 
or they may be points of higher order contact discontinuities. When compared with methods proposed by 
other investigators, the LSPE method was found to be at least as accurate and, often, significantly more 
accurate than other methods. Also, the LSPE method was found, in general, to be more robust and less 
sensitive to the effects of “closely spaced” singularities than other methods. 

However, there are some issues connected with the LSPE method that we feel need further study. For 
example, a good rationale for the choice of the parameter R in the definition of the objective function 

N 

( 10 . 1 ) E= J2 

j=N + l-R 

needs to be established. We have required that TV — R 1, in order for the asymptotic form of the Fourier 
coefficients to be valid, and that 2 R is greater than the number of parameters being estimated. Intuitively, 
the freedom to choose “larger” values of R has several advantages, including a smaller sensitivity of the 
LSPE method to “small errors” in the Fourier coefficients. In most of the results reported in this paper, 
the choice of R has been relatively ad-hoc, but some effort was made to keep it unchanged for each example 
considered. However, we observed that, in some cases (in Example 2, for instance), for different values of 
TV the convergence rate of the optimization routine seemed to be somewhat dependent on the choice of R . 
This sensitivity to the choice of R has yet to be studied in any detail and all of our studies in this regard 
are preliminary in nature. It does seem reasonable, though, to conjecture that this sensitivity is due to the 
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decreasing magnitude of the Fourier coefficients with increasing N. For large N and Af, the limitation of 
having only a finite amount of precision in the coefficients causes “large” relative errors in the coefficients 
that are used. On the other hand, the idea of using larger values of R , for larger N and M, in order to 
average out these relative errors, may not always be practical, as larger values of R may cause a slowing in 
the convergence rate of the optimization method. 

Also, a good criterion for the choice of the weights {wj} in the definition of E needs to be determined. 
In particular, the proper assignment of the weights {wj } is crucial in the design of an optimal least-squares 
optimization method. All the results reported using the LSPE method (unless otherwise noted) were obtained 
using Wj = j. Some experimentation with assigning the weights shows that this choice is not uniformly 
optimal However, all of our experimentation with {iUj} has been preliminary, and more study is required 
to determine the actual sensitivity of the LSPE method to the choice of Wj. 

Two important kinds of singularities in a function that are often encountered and which we have not 
addressed are algebraic and logarithmic singularities. Unlike Lipschitz functions, the Fourier coefficients of 
functions with these types of singularities do not have the same asymptotic form as for Lipschitz functions. 
Consequently, a new method to determine the locations and characteristics of the singularities of such a 
function need to be developed. Even after the singularities have been characterized, the basis functions 
(5 n (a;)} are no longer applicable, and some “new” basis functions with the appropriate algebraic or loga- 
rithmic singularity must be constructed. For an algebraic singularity, one possible candidate for a new basis 
function is 

oo 

G{x,p) = |sin(x/2)| p = y + cos 0' x )> 

3 = 1 

where, for p > — 1 , 

_ 2 r(( P + 1)/2) ( -i) j r(p+i) 

0 v/tF r(i +p/2)’ J 2* r(i+j+p/2) r(i-j+p/2)’ 3 - ' 

Here T ( z ) is the usual Gamma function. Near x = 0, we have 

c(*.p) = ! fr {i+o (**)}, 

and, hence, G(x,p) could be used as a “basis function” to simulate a p th order algebraic singularity in a 
function /. For a function with a logarithmic singularity, the function 

oo 

L(x) = - log |2 sin(x/2)| = V - cos (Jx), 

j=l J 

might serve as an appropriate basis function, since, near x = 0, 

L(x) = - log \x\ + 0{x 2 ). 

We also feel that many of the ideas we have presented for Fourier series can be extended to other 
orthogonal bases, as well. In particular, it is natural to expect “Gibbs-like” properties in “all” series using 
orthogonal sequences that approximate discontinuous functions. Gray and Pinsky [12] report a Gibbs like 
phenomenon for a Fourier- Bessel series of a piecewise smooth function, which displays a strange oscillation at 
the origin, quite unrelated to the local behavior of / at that point. We feel that the problem of constructing 
a high order accurate sequence of approximations using information from a finite set of general orthogonal 
series coefficients can be addressed using extensions of several of the methods we have presented. 


21 



In addition, we feel that the ideas presented here will find useful applications in several different areas. 
Some particular applications currently being pursued include numerical shock capturing, image resolution 
enhancement, and time series analysis. 
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Appendix A. Convergence of initial estimates. In this appendix we outline a proof of the rates 
of convergence as N — > oo of the estimates {x s } of the locations {#*} and the estimates of the 

magnitudes {[/(x 5 )]} of the simple discontinuities of /, obtained from the extrema of the partial sums , 

defined in Eq.(2.2). 

Using Eq.(2.2) we can write 

(A.l) Q Sm ^ du^j D n (x) = E sin (]^i) { bj cos U x ) ~ a j sia{jx)} . 


Then the condition that dD^/dx = 0 implies 

N 

(A.2) 

i=i 


0=EJ («i - (jx) + bj sin (j x)} . 

Using the asymptotic form (Eq.(3.1)) of the coefficients Eq.(A.2) can be written as 

n 

° = 51 [/(*.)] 


(A.3) 

where 


E sin {nTi) _ 

E [/'(*»)] jEjsm (j^lj cos 0 (*-a:.)) | +T n (x), 
t n (x) = 5Z sin ( i J {X cos U x ) + h sin 0 X )} . 

3 = 1 ' ' 

with dj = 0(j~ 2 ) and bj = 0(j -2 ), as j — > oo. 

Let Xfc = Xk + tfc denote the location of the extremum of Dn nearest to Xk- We then write Eq.(A.3), 
evaluated at x = x^, as 

° = [/ (^fc)i E sin ( J sin 0‘ e *) 

+E if (*.)] | E sin (atX“t) sin 0' A fc.») 
b =1 

- [/' (**)] E ) sin (atxt) cos (J e *=) 

“E^ {Ej™ (jv+j) -**■)} + T w(ifc), 


(A.4) 


where A^ s = + e* — for s ^ k. Since / is assumed to have only a finite number of points of 

discontinuity in the interval (— 7r,7r], we can assume that |A* |S | > A > 0. Then, using elementary methods, 
we can show that, for any A/0, 

N 


(A.5) 


^^ 8in (]vTi) sin(jA) = 0(Ar 1 


), m > 0, 
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(A-6) £ i sin ( jy^i) cos 0' A ) = 0(JV 1 ), m > 1, 

as AT — ■> oo. Using these results, we see that the second and fourth terms on the right side of Eq.(A.4) are 
each 0(N _1 ), as A -> cx). Also, T^(xk) = 0(N~ x ), as N — > oo, since the coefficients {dj} and are 
each 0(j ~ 2 ), as j — ► oo. If we now assume that jtk = o(l), as AT — > oo, for all j < AT, we can write Eq.(A.4) 
as 


(A.7) 

We now use the facts that 


0 = e k - (l + 0((iVe fc ) 2 )) ' [/(^fc)lX^Jsin - 

i=i ^ ^ 

[/' (**)] (1 + 0 ((iVe fc ) 2 )) • f; ^ Sin (^) + 

7 = 1 ^ V / 


(A.8) 


J = 1 7 


as AT — > oo, to write Eq.(A.7) as 


(A.9) 


- i r . iziM , m ^ r -3^ 


“ = (jv+tfX ^W +0(w_3) ' 


as TV — ► oo. Prom Eq.(A.9) and the definition of e*;, we see that x* = + 0(AT 2 ), as TV — ► oo, as asserted 

in Theorem 2.1. 

Also, using some of the same steps that led to Eq.(A.4), we can write 


(A.10) 


Si(fl-) ■ D N (x k ) = [/(a*)] 

+E [/(*.)] 

5=1 

s^k 

+£ [/'(*«)] 


| Ylj 1 sin ( j^r) | ( x + ° ((Ac*:) 2 )) 


5 — 1 


f N 

Ei 1 sin 1 
o =1 

( 3* N 

) cos(jA fc , a )| 

\N + 1, 

{gH 

' J* > 

| sin(jA fcjlS )| 

k N + l) 


Using Eqs.(A.5) and (A.6), we find that each of the terms beyond the first on the right side of Eq.(A.lO) is 
0(N ~ : ), as N — ► oo. Then, using Eq.(A.8), we find that Eq.(A.lO) can be written as 


D N (x k ) = [/(x fc )] + o(Ar- 1 ), 


as stated in Theorem 2.1. 
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Appendix B. Exponential convergence. We now outline a proof that the sequence 
converges exponentially to / in the maximum norm, outside the union of a finite number of small open 
intervals that contain the points of singularity of /. 

To establish this exponential convergence, we first define the error terms 


Em,n{x) = /(*) - = E$ n (x) + E™ n {x), 


?( 2) 


where 

(B.l) 


= f( x ) - EmIn( X ) = /m,/v(x) — Jm,n( x )' 


( 2 ) 


Here f\t,N is defined by the right side of Eq.(9.4) with x* and J^ s replaced by x s and [/^(x 5 )] , respectively, 
in the definitions (9.2)-(9.5). From the results of [6], E^\ M (x) decays to zero exponentially, as M — ► oo. 


(2) 

To show that E K M ; XM (x) also decays to zero exponentially, we note first that, from its definition, 

(B.2) 
where 

oo 

(B.3) Yk,ff(x) = E a *,j cos(jx) + bk,j sin(jx). 

j=N + 1 

Thus, it follows from Eqs.(B.l), (B.2), and (9. 2)- (9. 5) that 

M n 

e { m\ n (*) = EE Y ^ x - x «) - - <)) ■ 

k—0 5=1 

Here the constants are defined by Eq.(9.3) with a replaced by [/( fc )(x s )] . To facilitate our proof 

below, we rewrite E$ N (x) as 

(B-4) E%\ n {x) = E^{x) + E%% 

where 


N M n 

Im,n( x ) = + E a,j cos (jx) + bj sin (jx ) +EE- 4 *.* ~ x *) > 

j= i 


k=0 5=1 


M n 


(B.5) 

and 

(B.6) 


E m'n(. x ) = E E Ak ' s _ X >) ~ Y k,N(x ~ X*)} , 


k= 0 s=l 


*£$(*) = E E i A ^> - 4U Y ^ x - *:)■ 


M n 


fc=0 5—1 


In order to show that E { MXM (x) decays to zero exponentially, as M — ► oo, it suffices to show that both 
(2 11 (2 21 

( x ) \m( x ) decay to zero exponentially, as M — ► oo. 

Using Eqs.(6.3)-(6.4), we can write 


(B.7) 


(B.8) 


X B X s 




(14-0 (l/N )) , 


J* kl , - [f W (*.)] = (! + O (W) , as N - oo, 
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for s = 1,2, n, and k = 0, 1, M. Here an d are certain constants. 

We first outline a proof that E decays exponentially to zero. From Eqs.(9.1), it follows that 

(B-9) M < and |6 fc)i | < jgy , 

where C\ and C 2 are constants independent of k and j. Therefore, it follows from Eq.(B.3) that there exists 
a constant Ci, independent of k and N, such that 


(B.10) 


max 

— 7T<X<7T 


OO C \ 

K.'.ISft E foti>l, 


j=N+l • 


which follows easily by the integral comparison test. Using Eq.(B.3), Eq.(B.7), the integral comparison test, 
and Taylor’s theorem, we obtain 


(B.u) IW* - *:> - n.»(x - x,)| 

for s = 1,2, ..., n, k = 2, M, and for 


(B. 12 ) ja: — a:* | > 0(tp Ms N~ M ~ 2 ). 

Here, ^m,s = (1 + 0(1/AT)), and C 2 is a constant independent of M and k. We now assume that there 

exist some positive constants 8 1 , 82 , £ 3 , Ai, A 2 , and A 3 , such that the following bounds hold: 

\&M,k,s\ < Ai M 81 A:!, 


\Ak,s\ < A 2 8 % A:!, and 


(B.13) |^ M>S | < A 3 <5" Ml, for 3 = 1, ..., n. 

The restrictions (B.13) are mild, and are motivated by a similar study in [6]. Then, using Eqs.(B.ll) and 
(B.13), we find that, for A: > 2, 


cM \ri 

• in, *(*-*:)! <di$ n, 3 • 


N M + 2 (k- 1)’ 

where d\ is a constant. Using Stirling’s approximation for Ml and setting N — AM, we find 

M n /, , \M 

(B.i4) VV|^, s |.|n, JV (x-x 8 )-n, N (x-x;)i<(i 2 VM(^p) , 

fc=o*=i \ Ae / 


where d 2 is a constant independent of M. Thus, for any fixed values of J 2 and 83 (with <$ 2^3 > e), we can 
select A large enough so that (<5 2 ^3)/(^ e) < 1. Consequently, E^ l ^ M {x) decays to zero exponentially, as 
M — ► 00 . 

We now consider E^ 2 ^ M (x). From the definitions of A £ s and Ak yS , and Eqs.(B.8), it follows that 


(B.15) 


Als 




< 


k 

IgA fjcy* | \ a M,k-i,s\ 

jyM+l-fc + 2 ^ ^yM+l+i-fc 
i—2 
i even 


IA 


k,k-i\ 
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where aM, -i,s = 0 , /3 kk _ 2 = [5^ 2 (0)], and, in general, 0 kyk _ { is a function of the jumps (sj^(O)], j = 1, 
We now conjecture that the following bound holds for fi k fc „ i : 


<B16 > ^7)t(|)‘ 

Although we do not have an analytical proof for this bound, it has been verified using Mathematica [ 22 ] 
for 0 < i < k < 300. Therefore, using Eqs.(B.13) and (B.16), and assuming N > M/(e<5i), it follows from 
Eq.(B.15) that 


(B.17) 


^fc,s — ^k,a < 03 


M8\k\ 

]yM+\-k ’ 


where C 3 is a constant independent of M and k. Using Eqs.(B.17), (B. 10 ), and (B.13) we obtain 

i A . * | . . ... 7 M 5\ (k - 1 )! 

| A k % s - Afc,s| ■ \ Y k,N(x-X a )\ < d 2 J fAf+i * 

Therefore, using Stirling’s approximation for M! and setting N — AM, we find 


(B.18) 


M n / A \ 

£ E \ A *’» - A U • in,N(* - *:)\ < ( f- 

fc= 0s =i v Ae / 


M 


As a consequence, for a fixed <5 1 , if we select A such that A > <?i/e, E^ 2 ^ M (x) decays to zero exponentially, 
as M — ► oo. 

It now follows from Eqs.(B.4)-(B.6), (B.14) and (B.18) that 


(B.19) 


r (2) 
J M,\M 


(*) 




M 


where d\, d 2 are certain constants. From equation (B.19) we see that, for any fixed values of 8\, S 2 and <5 3 , 
if we select A so that 


(B.20) A>max(<5i/e, ^ 2 ^ 3 /e) , 

/Q\ 

then E k m \m(x) decays to zero exponentially, as M — ► 00 . 

Consequently, Em,\m{x) decays to zero exponentially in the maximum norm for x in the domain D. 
Here D consists of the interval [ — 7r, 7 t], with the union of a finite number of “small” open intervals, J a , for 
s = 1,2, n, surrounding the points of singularity of /, removed. In order to show that the total measure of 
the union of these “small” open intervals, 7 S , decays to zero exponentially, we observe from Eq.(B.12) that 
the length L s of the interval I a satisfies 


Ls < N - M - 2 ). 


Using the bound on V>m,s from Eqs.(B,13) we find that, for N = AM and for large M, 


L s <0 



Thus, it follows that the total measure, Y^=i °f the union of the intervals goes to zero exponentially, as 
M — ► 00 . 
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Finally, we show that yflf XM 0*0 j converges exponentially to f'(x) in the maximum norm in the domain 
D, as M — ► oo. We observe that, from the results of [ 6 ] and from Eqs.(B.l), (B. 4 ), (B. 5 ) and (B. 6 ), it suffices 
to show that both ( x )) an d decay to zero exponentially, as M — > oo. Using Eqs.(B. 3 ), 

(B.13), the integral comparison test, and Taylor’s theorem it follows that 


(B.21) 

and 




< B22 > 

for k = 3, M. Here, C 4 and C 5 are constants independent of k and M. Therefore, using Stirling’s approx- 
imation for Ml and setting N = AM, it follows that 


(B.23) 

and 

(B.24) 






(«„(*))'| < M 5 / 2 


M 


Thus, assuming that A satisfies the condition (B. 20 ), it follows that aa/M) converges to f f (x) expo- 
nentially in the maximum norm in the domain D, as M — ► 00 . Using mathematical induction, it follows 
that the first l derivatives of converge exponentially to the corresponding derivatives of /(x), as 

M — > 00 . This completes our proof. 
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Fig. 1. Plots of Djy(x) for Example 1 using N = 16 (dashed line) and N = 64 ( solid line). 
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Fig. 2. The normalized errors N 2 \x\ - Xi\ (solid line), obtained from D N (x), and N 2 |xj -xi| (dashed line), obtained 
from the LSPE method , for Example 1, plotted as functions of 1/N. 


31 








Fig. 4. A surface plot of the error function E (see Eq.(3.3)) for Example 1, as a function of X\ and Jo,i } near the 
minimum of E. Here N = 32, R— 12, and u>j = j. 
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Fig. 5. Plots of D i,n(x) for Example 1 using N = 32 (dashed line) and N = 64 (so/td fine/ 
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Fig. 6. The reconstructed approximation /* 32 (x) (dashed line) for the function f(x) (solid line ) of Example 1, using the 
parameter estimates obtained by the LSPE method with M ~ 2 and N = 32. 
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Fig. 8. The reconstructed approximation (/ 2 * 32 (i))" ( dashed line) for the function f"(x) (solid line) of Example 1, 
the parameter estimates obtained by the LSPE method with M — 2 and N = 32. 
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